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Recent progress in simulation methodologies and in computer power allow first principle 
simulations of condensed systems with Born-Oppenheimer electronic energies obtained by 
Quantum Monte Carlo methods. Computing free energies and therefore getting a quantitative 
determination of phase diagrams is one step more demanding in terms of computer resources. 
In this paper we derive a general relation to compute the free energy of an ab-initio model with 
Reptation Quantum Monte Carlo (RQMC) energies from the knowledge of the free energy of 
the same ab-initio model in which the electronic energies are computed by the less demanding 
but less accurate Variational Monte Carlo (VMC) method. Moreover we devise a procedure 
to correct transition lines based on the use of the new relation. In order to illustrate the 
procedure, we consider the liquid-liquid phase transition in hydrogen, a first order transition 
between a lower pressure, molecular and insulating phase and a higher pressure, partially 
dissociated and conducting phase. We provide new results along the T = 600/^ isotherm 
across the phase transition and find good agreement between the transition pressure and 
specific volumes at coexistence for the model with RQMC accuracy between the prediction 
of our procedure and the values that can be directly inferred from the observed plateau 
in the pressure-volume curve along the isotherm. This work paves the way for future use 
of VMC in first principle simulations of high pressure hydrogen, an essential simplification 
when considering larger system sizes or quantum proton effects by Path Integral Monte Carlo 
methods. 

Keywords: Monte Carlo Methods; Quantum Monte Carlo; High pressure hydrogen. Free 
energies and phase diagrams 



1. Introduction 

First-principles simulation has become an essential method to investigate the phys- 
ical behavior of condensed matter systems, in particular when the chemical nature 
of the interaction among the nuclear degrees of freedom depends on the external 
conditions, like for instance in systems under extreme conditions of pressure and 
temperature. Those methods generally treat systems of classical nuclei at finite 
temperature and electrons in their instantaneous ground state in the framework of 
the Born-Oppenheimer (BO) approximation, although extensions to treat quantum 
nuclei at finite temperature within the Path Integral formalism, and/or electron 
at finite temperature with fractional (fermi) occupation have been developed. The 
electronic structure at fixed nuclear positions is generally computed within the 
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Density Functional Theory (DFT) framework [l| or, when higher accuracy is re- 
quired, by Quantum Monte Carlo methods (QMC) DFT is, in principle, an 
exact theory but, in practice, it is based on uncontrolled approximations for the 
exchange-correlation functional. Despite some well documented limitations, DFT 
is often accurate and fast enough to be used in conjunction with Molecular Dy- 
namics to perform dynamical studies for systems of several hundred nuclei and 
extract physical information (first-principles Molecular Dynamics, FPMD) 
Ground state QMC comes in two different flavors: Variational Monte Carlo (VMC) 
and Projection Monte Carlo (either Diffusion Monte Carlo (DMC) or Reptation 
Quantum Monte Carlo (RQMC)). VMC is faster but the accuracy of its predictions 
is limited by the skill in designing and optimizing suitable many-body trial wave 
functions. RQMC is a method to automatically optimize trial wave functions; it 
provides the most accurate description of the electronic properties but has larger 
computational requirements (about an order of magnitude larger than VMC, not 
counting the VMC optimization process). It is an exact method for bosons since 
it is able to provide the properties of the bosonic ground state in a time which 
is polynomial in the number of particles. In the case of fermions (electrons) the 
scaling is exponential and for practical calculations we have to resort to the "fixed 
node approximation". This approximation makes the method variational with re- 
spect to the trial nodes; usually the results remain highly accurate if good trial 
nodal surfaces are employed 0]. With recent progress in optimization procedures 
0] the location of the trial nodal surfaces, suitably parametrized, could also be 
optimized to achieve higher accuracy. 

In recent years we have developed an ab-initio method in which electronic ener- 
gies are computed by ground state Quantum Monte Carlo (QMC) methods while 
nuclear degrees of freedoms are sampled by Metropolis Monte Carlo. This method, 
called Coupled Electron- Ion Monte Carlo (CEIMC) 0], has been successfully ap- 
plied to high pressure hydrogen. The results are a benchmark of DFT calculations 
in the region of phase diagram of interest for planetary physics [1, 0] and across 
the metal-insulator transition region in the fluid phase [id. 111]. 

The general picture emerging is that FPMD for high pressure hydrogen is gen- 
erally accurate away from the metal-insulator transition while QMC accuracy is 
required when approaching the metallization of the system. 

In CEIMC we sample the equilibrium distribution for the protonic degrees of 
freedom by a Metropolis based Monte Carlo scheme employing either VMC or 
RQMC to compute the difference between the BO energies of two nearby protonic 
configurations in the Metropolis method. The electronic calculation needs to be 
repeated for each attempted move of the protons and it represents by far the 
largest part of the computational load of the method. It is therefore very tempting 
to exploit VMC rather than RQMC in CEIMC. On the other hand if the trial wave 
function is not accurate enough, the use of VMC could provide biased results which 
can lead to inaccurate predictions of the physical behavior of the system [l3| ■ 

The key quantity in tracing transition lines is the free energy difference (either 
Helmholtz or Gibbs) between coexisting phases which is usually obtained by com- 
puting the absolute free energies of the coexisting phases. This can be achieved by 
the Coupling Constant Integration (CCI) method, a general strategy introduced by 
Kirkwood to adiabatically transform the system into another system of known 
free energy. This strategy is often used in ab-initio simulations to compute the free 
energy of the ab-initio model starting from the knowledge of the free energy of an 
effective classical representation of the same system, for which standard free energy 
methods can be easily applied [HI, 1^ lit. 

In this paper, within the framework of the CCI and of the CEIMC method. 
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we derive a relation which ahows us to obtain the free energy of the system with 
RQMC electronic energies from the knowledge of the free energy of the system with 
VMC energies. The simple relation, derived in the next section, is of practical rele- 
vance when tracing transition lines since it allows us to obtain the transition lines 
with RQMC accuracy by performing free energy calculations with VMC electronic 
energies. This strategy improves the efficiency of the method by roughly 10 — 20 
times, depending on the accuracy of the trial wave function. 

To illustrate the use of the new relation we consider the liquid-liquid phase 
transition in high pressure hydrogen, predicted earlier using the so-called 

chemical model, but only recently put on firm basis by first-principle simulations 
[nl, [H]. This is a first order transition between an insulating, molecular liquid 



phase at lower pressure and a metallic, partially dissociated liquid phase at higher 
pressure. It is predicted to occur below 2000K; on increasing pressure the specific 
volume is discontinuous, corresponding to a partial dissociation of molecules and 
to an abrupt jump of roughly four orders of magnitude of the DC conductivity [lH . 
Above T=2000K the process is observed to become continuous and a critical point 
is predicted to be located at = 2000K and Pc{DFT) = 85GPa, PdQMC) = 
120GPa. The transition line T(P) has a negative slope, meaning that lowering the 
temperature requires higher pressure to cross the transition, and it is predicted to 
meet the melting line of the molecular insulating crystalline structure (Phase I) at 
a triple point located at Tt = 700K and Pt = 220GPa within the DFT framework 
and Tt = 550K and Pt = 290GPa using CEIMC [lH. The latter prediction was 
based on an extrapolation to high pressure of lower pressure results. 

In this paper we performed new CEIMC calculations of high pressure hydrogen 
for T=600K and across the liquid-liquid phase transition. We performed calcula- 
tions using both VMC and RQMC with the aim of testing the accuracy of VMC 
across this particular phase transition. Having at our disposal the equation of state 
obtained by the two methods we can test the reliability of our new relation for 
correcting free energies. 

The paper is organized as follows. In the next section we will set the formal 
framework of RQMC and we will derive the formula to correct the free energy. In 
section [3l as an application of our strategy, we will present recent results obtained 
for the liquid-liquid transition line of high pressure hydrogen. Finally in section U 
we will draw some conclusions. 



2. Theoretical framework 

Let us consider a system of N monovalent ions of mass M and N electrons of 
mass m in a volume V and at thermal equilibrium with a bath at temperature T. 
If we assume the validity of the Born-Oppenheimer approximation and consider 
electrons in their ground state, the electrons provide a many-body potential for 
the ionic motion. The ionic "Hamiltonian" is 

n.m,^.}) = E ^ n ^ ^ ^"""^^ 

1=1 i-^j 

where {Pi} and {Rj} are the set of momentum and coordinate operators of the ions, 
respectively, and E{{'Ri}) is the electronic energy which depends parametrically 
on the nuclear positions. The thermodynamics of the system is derived from the 
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knowledge of the Helmholtz free energy F{T, p) 

F{T, p) = -keT log Wr (e-^^' ) (2) 



where /3 = l/ksT and the trace is over nuclear degrees of freedom. The term 
£'({Ri}) can be computed at various levels of approximation, for instance us- 
ing empirical effective potentials or Density Function Theory (DFT) or Quantum 
Monte Carlo (QMC) methods. Within the latter class of methods we can distin- 
guish between Variational Monte Carlo (VMC) and Reptation Quantum Monte 
Carlo (RQMC) which provide different level of accuracy. 

In VMC we make an ansatz for the many-body wave function ^^-(r^) = {r^\^|)T) 
and compute the electronic energy as a statistical average of the local energy 
Eioc{r^) = [ne4'T{r^)]/^T{r^) over the distribution |?/'T(r^)lV / IV-TCr^)!^ 

^ {iPt\iPt) /|VT(r^)p(ir^ / |VT(r^)Pdr^ ^ ' 

where = (ri, ■ ■ ■ ,rAr) is the vector of the 3N-electronic positions and Tie indi- 
cates the electronic hamiltonian 

In an electron-ion system, although not indicated explicitly, the trial wave function 
depends parametrically on the ionic coordinates and generally on parameters that 
need to be optimized by the use of the variational theorem. 

In RQMC the trial wave function is automatically optimized by the application 
of the projection operator e~*^= which filters out its excited states components. All 
excited states will be suppressed exponentially fast with increasing the rate of 
convergence to the ground state depending on the energy gap between the ground 
and the first excited state non-orthogonal to the trial function. Here t is a parameter 
with the units of inverse energy; we will refer to it as the projection time because 
of the analogy between the projection operator e~*^= and the propagator of the 
real time dynamics of the system, g"**^"/^. The total energy function of time is 
defined as 

mmmm) (VTie-^^Hv-T) (Vrie-^^nv'T) - ° 

n 

where i?o is the exact ground state □. 

Similar to a thermal partition function, let us define the generating function of 
the moments of "He at time t as 

Z(t) = (tAtIc-^^^I^t). (6) 
The total energy at time t is simply the derivative of the logarithm of 



^In case of fermions one must employ the fixed node approximation to avoid the sign problem and Eq is 
the fixed node ground state energy which is an upper bound of the true ground state energy 
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and the variance of the energy is the second derivative 

aUt) = {ine-Eit)f) = -^^Eit)>0 (8) 

which is non-negative by definition and implies that the energy decreases mono- 
tonically with time. The (fixed node) ground state is reached at large time (much 
larger than the inverse gap) and 

hm E{t) = Eq (9) 
lim a|(i) = (10) 

The last relation is the generalization of the zero variance principle to RQMC. Note 
that the average in eq.® is over the electronic ground state, not to be confused 
with the trace over ionic degrees of freedom in eq. ([2]). 

The VMC level of description is obtained at t = 0. This observation allows us to 
derive an expression to connect the free energy of the system with VMC electronic 
energy to the free energy of the system with RQMC electronic energy. Let us call 
T-Lii^t) the nuclear hamiltonian of eq. ^ with electronic energy E({'Ri},t) obtained 
after a time t projection, and Ft{T,p) the corresponding free energy. The time 
derivative of the free energy is 



dFt _ fdE{{Ki},t) 



dt V dt 



Mt (11) 



where (• • • ) indicates the trace over ionic degrees of freedom with statistical weight 
g-PHi (t) ^ fpj^g £j,gg energy at any positive time t is obtained from the knowledge of 
the VMC free energy Fo{T,p) as 

Ft{T, p) = Fo{T, p) - f ds{al{s))s (12) 
Jo 

where the second term on the r.h.s. is the integral of the average, over ionic sam- 
pling, of the variance of the electronic energy. Since the variance is positive by def- 
inition, the time derivative of the free energy is negative, implying that the free en- 
ergy is a monotonically decreasing function of the projection time t: Ft < Eq, Vt > 0, 
that is the free energy of the VMC system is an upper bound to the fixed node 
free energy which is obtained at infinite projection time. 

Eq. (jl2p is the main result of our work. It is particularly useful to increase the 
level of accuracy of transition lines from VMC to RQMC. Indeed once a transi- 
tion point has been detected with the VMC energies, RQMC simulations at fixed 
thermodynamic points inside the two competing phases can be run for increasing 
projection time and the corrections to the free energy of both phases can be com- 
puted as we will show in the next section. This correction will determine how the 
transition point moves in going from VMC to RQMC accuracy. This method is 
much faster than performing thermodynamic integration directly within RQMC. 

Note that, being a generalization of the CCI method, eq. (jl2|) holds only if the 
system does not cross a phase boundary during the time projection. 
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3. Results for high pressure hydrogen 



To test the use of eq. p2|) , we consider the first order Uquid-hquid phase transition 
in high pressure hydrogen. We have extended previous CEIMC calculations [ll|] to 
the T = 600K isotherm, in the range of pressures P € [220, 330] GPa correspond- 
ing to a range of specific volumes v € [9.42, 11.49]aQ and of = (3u/47r)^/^ G 
[1.31, 1.40] where we expect the transition to be located. We performed CEIMC 
simulations in the NVT ensemble, on systems of 54 classical protons and 54 elec- 
trons. We adopted backflow-Slater-Jastrow electronic trial wave functions j20l. l2ll]. 
with analytical RPA and numerically optimized two-body correlations and back- 
flow terms. The single particle orbitals of the Slater determinants were obtained by 
self-consistent DFT-LDA calculations. To reduce the flnite size effects we adopted 
Twist- Aver aged Boundary Conditions [12] with a flxed grid of 64 twist angles. 
Moreover, finite size corrections were added to energy and pressure 23, [13]. More 



details can be found in ref. [ill] 



On the left panel of Fig. [T] we show the pressure versus the specific volume along 
the T = 600-ftr isotherm, as obtained by CEIMC with VMC energies (blue squares) 
and RQMC energies (pink circles). In the RQMC calculations the total projection 
time t (t = 1.00H~^) and the time step t {t = 0.02ff~^) were chosen to pro- 
vide converged energies, while pressures were corrected for finite t and nonzero 
r. Numerical values are summarized in Table [H The two curves have the same 
qualitative behavior and are very close especially at large specific volume, where 
pressures from the two methods agree within statistical errors. Both curves have 
a plateau at P ~ 280GPa, signifying a weakly first order phase transition (diver- 
gence of the isothermal compressibility). A peculiarity of this transition, probably 
related to its weakly first order character, is the absence of metastable states which 
typically appear in other first order transitions, such as, for instance, the gas-liquid 
transition, and makes difficult to obtain a quantitative location of the transition 
point without computing the free energies of the competing phases. It follows that 
the coexistence pressure and the specific volume of the coexisting phases can be 
qualitatively located directly from the P{v) behavior without resorting to the free 
energy methods. 

The values of the coexistence pressure, inferred from the position of the plateau 
in the P{v) curves in Figure [H are of about ^(J^qmc) ~ 280GPa and P(vmc) ~ 
275.5GPa for the RQMC and the VMC data respectively. The EOS for the two 
different phases can be fitted with polynomial functions. At pressures below the 
transition, the system is in its normal phase (phase I), molecular and insulating. 
The curves P{v) can be well represented by a second order polynomial Pi{v) = 
PI + qiv + riv'^, with pi = 2604GPa, qi = -382.5GPa/a|] and r/ = 15.2GPa/a^ 
for RQMC data and pi = lOOSGPa, qi = -91. DGPa/aj] and rj = 2MGPa/af, 
for VMC data. For the phase II, the conducting fluid, a flrst order polynomial 
Pii{v) = Pi I + qiiv is sufficient to describe the pressure behavior; the values of 
the parameters are pu = 1143GPa and qn = — 86.2GPa/ag for RQMC data and 
pii = 1256GPa and qn = — 98.6GPa/ao for VMC data. The specific volumes 
of the two phases at coexistence can be estimated by extrapolating those fits up 
to the coexistence pressure, as shown in the right panel of Figure [TJ The volume 
change at transition is in both methods ~ 4%. The specific volumes of the phases 
at coexistence are ^^//(rqmc) ~ 10.00(4)aQ and ^^/(rqmc) ~ 10.31(4)aQ for RQMC 
data and ?^//(vmc) ~ 9.94(2)aQ and ■^/(vmc) ~ 10.30(2)aQ for VMC data. 

To test the ability of the correction formula (jl2p to predict the RQMC coexistence 
pressure of the metallization transition at T = 600i^r we proceed as follows. 
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a) We assume that the coexistence value for the pressure and the specific vol- 
umes from the qualitative visual analysis above are correct and provide our 
reference results for the VMC system. With this assumption and using the fits 
to the EOS's above, we reconstruct the Helmholtz free energy curves for the two 
phases up to a constant value, by integrating the VMC-EOS's. For a given phase 
a = I,II the free energy per particle at the specific volumes v is related to the 
free energy per particle at coexistence volume by 



fa{v) - faK) = - 




= -P.iv - <) - y - Kf) - y - (13) 

where the final expression has been obtained by integrating a second order poly- 
nomial fit of the equation of state. The free energies of two coexisting phases can 
be referred to a single reference value noticing that at coexistence 

fiiv}) - fiMi) = P^h - vj) (14) 

The resulting curves are plotted in Figure [21 The magenta open squares repre- 
sent the free energy per particle of phase II and the blue open circles the free 
energy per particle of phase I. In order to enhance their curvature and make 
the visualization of the common tangent easier, we subtracted the same linear 
common tangent £{v) = P^iv - vj) = 0.00934i//a|] v - 0.09303i7 from each 
curves. This common tangent (represented by the horizontal black line in Figure 
[2|) is, by construction, the coexistence pressure corresponding to the plateau in 
the VMC pressure curve as discussed above. 

b) Once the free energy curves for the VMC system are known we can apply the 
correction term in eq. (|12p to obtain the RQMC free energy curves. The most 
convenient strategy is to compute the correction at a single fixed specific volume 
for each phase Va and to reconstruct the RQMC-free energy curves by integrating 
the RQMC-EOS of the two phases. To this aim, at a single fixed density in each 
phase we perform several RQMC calculations for increasing projection time to 
numerically compute the correction term 

Saiva) = /i^«*'^(^.) - fr'^'iva) = lim /* ds{aUs,Va)), (15) 

The values of the reference volumes Va must be taken far enough from the coex- 
istence region to avoid crossing the phase boundary during the time integration 
process. To satisfy this requirement we have chosen vjj = 9.634aQ (r^ = 1.32) 
for the metallic liquid (phase II) and vj = 11.249aQ (r^ = 1.39) for the insulating 
liquid (phase I) . In Figure [3] we compare the proton-proton radial distribution 
functions g{r) at the two densities as obtained by CEIMC with VMC and RQMC 
electronic energies. We can see that at both densities the system is in the same 
phase regardless of the electronic method used. At vj = 11.249aQ both g{r) 
present a sharp peak centered at r ~ 1.4ao, which indicates that the system 
is in the molecular phase; at vjj = 9.634aQ the molecular peak has completely 
disappeared and the form of g{r) is typical of a monoatomic fluid. 

In order to evaluate the integral in eq. (|12p . at each reference specific volume 
we performed additional CEIMC-RQMC simulations with projection times t = 
0.08, 0.20, OAOH~^ and a time step r = 0.02H~^. The behavior of a^/N with t, 
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corrected by finite time step error, is sliown in Figure U] and is well fitted by a 
stretched exponential function: s{t) = ae"^*/''^^ and its time integral is obtained 
by using the relation 

f .-<./'r.*.^r(i) (16) 

where T{x) is the Gamma function. The values of the fitting parameters in 
the two phases are a/ = 0.01422if^ bj = 0.0440if-\ cj = 0.797 at v = vi and 
an = 0.01042ij2^ bn = 0M55H-\ cn = 0.987 at t; = vn and the RQMC-VMC 
Helmholtz free energy difference per particle results to be dj{vi) = —0.71{l)mH 
for the molecular insulating phase (phase I) and 6n{vii) = — 0.48(l)mi7 for the 
metallic phase (phase II). 

c) finally we use the RQMC-EOS to construct the RQMC free energy curves 
and to determine the coexistence pressure and volumes between the phases I 
and II by the common tangent construction. In figure [2] the RQMC free energy 
curves of phase I and II are represented by the closed turquoise circles and the 
closed purple squares, respectively. Despite the fact that the correction from the 
VMC free energy is quite smal0, less than ImH per particle over the range of 
volumes under analysis, it has a significant effect on the coexistence: the common 
tangent construction applied to those curves provides a coexistence pressure of 
P'^ = 281(l)GPa and specific volumes of the two phases at coexistence of v'f = 
10.32(4)a|] and vj^ = 9.99(4)a|]. Those values are in excellent agreement with 
the values directly inferred from the RQMC pressure curve {P'^ = 280(l)GPa, 
= 10.31(4)aQ and v^^ = 10.00(4)aQ), being the differences between the two 
predictions within error bars. 

4. Discussion and Conclusion 

Within the Coupling Constant Integration strategy and for the Coupled Electron- 
Ion Monte Carlo method, we have derived an original relation to obtain the free 
energy of an ab-initio system with Reptation Quantum Monte Carlo electronic en- 
ergies from the knowledge of the free energy of the same system but with Variation 
Monte Carlo electronic energies, a considerably less demanding calculation. This 
relation will be useful in improving the accuracy of the CEIMC method in tracing 
transition lines, therefore predicting phase diagrams, without increasing too much 
the already large computational requirements of the method. In order to test the 
use of the new relation and to illustrate the procedure to correct transition points 
we have considered the liquid-liquid transition in high pressure hydrogen, a first 
order transition between a lower pressure, molecular and insulating phase (I) and 
an higher pressure, partially dissociated and conducting phase (II). We have re- 
ported new CEIMC results along the isotherm at T = GOOi^T in the density range 
corresponding to € [1-31, 1.40] across the transition, an isotherm not considered 
before. Not presenting metastable states, this transition can be detected by the 
appearance of a plateau in the pressure-volume curve and the values for coexis- 
tence pressure and for the specific volumes of the phases at coexistence can be 
inferred without performing free energy calculations. Nonetheless, as detailed in 



^We can have an idea of the relative amount of the free energy corrections by comparing the difference 
^/('^i) ~ = 0.23(2)m/f with the VMC free energy difference between phase f and II at the same 

thermodynamic points, fi{vf) — fjiivn) ~ li.bmH: the contribution of the correction is then quite small 
as expected, about 1.5% of the VMC free energy differences. 
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the paper, it can still be used to illustrate our procedure for correcting transition 
lines. In particular, we get excellent agreement between the RQMC critical pressure 
inferred from the data {P'^ = 280(l)GPa) and the value obtained by our procedure 
[P'^ = 281(l)GPo). A similar agreement is obtained for the specific volumes of the 
two phases at coexistence. 

The close agreement between VMC and RQMC based Equation of States (EOSs) 
and transition points results from the high quality of our trial wave function for high 
pressure hydrogen even across the metal-insulator transition, a notoriously difficult 
region for DFT based first principle methods. This agreement strongly supports 
the use of VMC in CEIMC investigations of the phase diagram of hydrogen at high 
pressure, allowing calculations on larger systems or systems with quantum protons 
within the Path Integral formalism |3] . Moreover the formula derived in the present 
work can be used to improve the CEIMC- VMC transition lines to RQMC accuracy. 
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Tabic 1. Pressures along the T — 600/^ isotherm, for RQMC (third column) and VMC (fourth column) data 
sets. In the first column, the values of rg and in the second the corresponding specific volumes v are reported. 



rs 




PRQMc{GPa) 


PvMc{GPa) 


1.310 


9.417 


331.6(6) 


327.2(6) 


1.320 


9.634 


311.4(4) 


305.2(6) 


1.330 


9.855 


293.9(6) 


284.0(6) 


1.340 


10.08 


283.8(7) 


275.9(6) 


1.350 


10.31 


277.2(7) 


275.0(6) 


1.355 


10.42 


277.8(8) 


- 


1.360 


10.54 


271.5(6) 


262.8(6) 


1.370 


10.77 


253.1(6) 


252.6(6) 


1.380 


11.01 


240.7(6) 


239.2(5) 


1.390 


11.25 


231.0(6) 


229.6(6) 


1.400 


11.49 


220.4(6) 


217.8(6) 



1.31 1.33 1.35 1.37 1.39 1.41 



1.31 1.33 1.35 1.37 1.39 1.41 



340 



320 



300 - 



280 



D- 260 - 



240 



220 - 



200 




Q. 




1.31 1.33 1.35 1.37 1.39 1.41 





340 






320 






300 




n3 






CL 


280 




C5 


260 




Q. 


240 






220 






200 






Figure 1. Pressure P versus specific volume v, along tlie isotlierm T = 600K for VMC and RQMC data. 
Statistical error bars are smaller than the symbol size. Left panel: RQMC (pink circles) vs VMC (blue 
squares) comparison. Lines arc guides to the eyes. Right panels: polynomial fits (dashed black lines) of 
the equations of state for RQMC (top panel) and VMC (bottom panel) data. The estimated coexistence 
pressure is represented by the horizontal dotted line, corresponding to P = 275.5GPa for the VMC curve 
and P = 280GPa for the RQMC curve. The volumes of the two phases at coexistence are estimated as the 
intersection of this line with the polynomial fits and are Vjj ~ 10.00(4)a|j and Vj ~ 10.31(4)aQ for RQMC 
data and v'jj ~ 9.94(2)ag and ~ 10.30(2)ag for VMC data. 
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Figure 2. Common tangent construction applied to the specific Helmholtz free energy curves for the system 
with VMC energies (open symbols) and to the corrected curves (closed symbols). To enhance the curvature 
a linear term, £{v) = 0.00934?; — 0.09303, has been subtracted to all curves. The common tangents are 
represented by the black lines. After the addition of the correction the coexistence pressure changes from 
P" = 275.5(5)GPa to P" = 281(l)GPa. The qualitative coexistence pressure inferred from the RQMC 
data is P" ~ 280(l)GPa. 





Figure 3. Comparison between the proton-proton radial distribution functions g{r) obtained with the two 
electronic methods at rs = 1.32 (metallic phase (II), loft panel) and = 1.39 (molecular insulating phase 
(I), right panel). 
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Figure 4. Average variance of the BO energy per particle a'^/N as a function of the total projection time 
t, at T = 600-fC and for two densities at the two edges of the liquid-liquid transition: Ts = 1.39 (blue circles, 
molecular insulating phase I) and = 1.32 (pink squares, metallic phase II). The behavior ofa^/N with 
t can be fitted as a stretched exponential s{t) = ae"''^" , with a = 0.01422//^, b = 12.05//, c = 0.797 at 
rs = 1.39 and a = 0.01042J/2^ b = 21.1H, c = 0.987 at = 1.32. Values at finite t are corrected for the 
nonzero time step error. 



